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We investigate thermodynamic properties of anisotropic colloidal dumbbells in the frameworks provided by the 
Reference Interaction Site Model (RISM) theory and an Optimized Perturbation Theory (OPT), this latter 
based on a fourth-order high-temperature perturbative expansion of the free energy, recently generalized to 
molecular fluids. Our model is constituted by two identical tangent hard spheres surrounded by square-well 
attractions with same widths and progressively different depths. Gas-liquid coexistence curves are obtained by 
predicting pressures, free energies, and chemical potentials. In comparison with previous simulation results, 
RISM and OPT agree in reproducing the progressive reduction of the gas-liquid phase separation as the 
anisotropy of the interaction potential becomes more pronounced; in particular, the RISM theory provides 
reasonable predictions for all coexistence curves, bar the strong anisotropy regime, whereas OPT performs 
generally less well. Both theories predict a linear dependence of the critical temperature on the interaction 
strength, reproducing in this way the mean-field behavior observed in simulations; the critical density - that 
drastically drops as the anisotropy increases — turns to be less accurate. Our results appear as a robust 
benchmark for further theoretical studies, in support to the simulation approach, of self-assembly in model 
colloidal systems. 


I. INTRODUCTION 

In recent years, a considerable number of experimen- 
tali“— and numerical^i^— studies have been devoted to the 
investigation of phase behavior and self-assembly proper¬ 
ties of colloidal dumbbells. Specific interest in this class 
stems essentially from the possibility to accurately tune 
the aspect ratio of the constituting spheres, as well as 
their interaction properties, so to obtain rich and fasci¬ 
nating phase behaviorsIn particular, if one of the 
two particles is solvophilic and the other one is solvo¬ 
phobic, colloidal dumbbells represent a simple exam¬ 
ple of surfactant (Janus dumbbells) constituting a 
molecular generalization of Janus spherical particles 
largely investigated in the last years because of their pe¬ 
culiar self-assembly properties and phase behaviorsi^^— 

While a large number of experimental and numerical 
investigations have involved colloidal dumbbells, consid¬ 
erable less attention has been paid to theoretical stud¬ 
ies: the framework of potential energy landscape has 
been adopted to describe the self-assembly of dumbbells 
into helicesj^ the competition between self-assembly and 
phase separation has been documented by integral equa¬ 
tion theories)^ the same approach^Si^ and fundamental 
measure theory^^ have been adopted to investigate struc¬ 
ture and thermodynamics of hard dumbbells. Notwith¬ 
standing all such studies, a systematic theoretical de- 
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scription of the phase behavior of colloidal dumbbells is 
still lacking. In order to fill this gap, here we calculate 
systematically free energy, pressure, chemical potential 
and phase diagram of anisotropic colloidal dumbbells, 
so to investigate how such properties change by tun¬ 
ing the interaction potential. We consider to this task 
a series of model dumbbells constituted by two identi¬ 
cal tangent hard spheres surrounded by square-well at¬ 
tractions with same widths and progressively different 
depths. We carry out our investigation in the frameworks 
provided by the Reference Interaction Site Model theory 
(RISM)^ and a recently proposed Optimized Perturba¬ 
tion Theory (OPT)4^i2^ Theoretical predictions for the 
gas-liquid phase coexistence are compared with previous 
Monte Carlo results by usi^ and other authors<^ The 
present work extends our preliminary OPT analysis of 
the symmetric case^ and our RISM study on a slightly 
different dumbbell modeli^ 

As for the RISM theory^^^ it constitutes a molecular 
generalization of the Ornstein-Zernike theory of simple 
fluids4i Originally devised for rigid molecules consti¬ 
tuted by hard spheres^ the theory was later extended 
to more complex systems, like colloidal models; in partic¬ 
ular, it was used to investigate the thermodynamic and 
structural properties of discotic lamellar colloids j^^do the 
self-assembly in diblock copolymers modeled as “ultra- 
soft” colloids;^ the interaction between colloidal parti¬ 
cles and macromoleculesthe crystallization and solva¬ 
tion properties of nanoparticles in aqueous solutions;^ 
the liquid structure of tetrahedral colloidal particles^, 
and the self-assembly properties of Janus rods4^ 
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FIG. 1. Schematic representation of colloidal dnmbbells in¬ 
vestigated in this work. Starting from the symmetric case 
£ii = £12 = £22 (top, left), the interaction becomes pro¬ 
gressively asymmetric by reducing £11 = £12 = 0.7, 0.5, 0.3, 
0.l£22. The weakening of £11 is signaled by progressively “fad¬ 
ing” the color of the first site. 


by two identical tangent hard spheres of diameter cr, la¬ 
beled as 1 and 2, interacting with sites 1 and 2 of an¬ 
other dumbbell at distance r via square-well (SW) at¬ 
tractions with same width (A = 0 . 5 ) and different depths 
£ij (h j — Ij 2). 

{ 00 if r < cr 

-£ij if a <r < a + Xa ( 1 ) 

0 otherwise. 

Parameters a and £22 provide, respectively, the unit of 
length and energy, in which terms we define the reduced 
temperature T* = /cbT/£ 22 (with fee as the Boltzmann 
constant), number density p* = pa^ and pressure P* = 
Po'^/e-2,1- In all calculations we have fixed £i2(= £21) = 
£11; then, we have studied a sequence of cases obtained by 
progressively reducing £11 starting from the symmetric 
case, £11 = £22, to £11 = 0 . 7 , 0 . 3 , 0.1 (in units of £22). 


As for the perturbation theory, this scheme played a 
significant role in describing the thermodynamic prop¬ 
erties of a large variety of fluids 4 S“— The theory was 
originally proposed as a high-temperature expansion by 
several groups and then recast in a very efficient com¬ 
putational tool by Barker and Hendersoni^^— It has 
been recently adapted to cope with Janus^ and tri-block 
Janus^ fluids. Here we use an optimized version of this 
theory based on a fourth-order high-temperature per¬ 
turbative expansion of the free energy, recently devel¬ 
oped to fit a large set of simulation data for square-well 
fluids.— Recently, one of the authors generalized this ap¬ 
proach to deal with more complex systems, like molecular 
fluids— Within this formalism, the effective energy and 
range of square-well interactions depend on the molec¬ 
ular anisotropy. In this way, it is possible to calculate 
analytical equations of state for a large variety of molec¬ 
ular models like for instance dimers, trimers, chains and 
spherocylinders— 

The importance of assessing the accuracy of these the¬ 
oretical approaches is crucial, as a large number of de¬ 
manding numerical calculationsi^ is typically required 
to investigate the experimentally relevant heterogeneous 
dumbbells casei^ A reliable theoretical tool would then 
help to focus on those systems of great interest. 

The paper is organized as follows: in Section II we 
provide the details of interaction potential and RISM and 
OPT schemes. Results are presented and discussed in 
Section III and conclusions follow in Section IV. 


II. MODELS AND THEORIES 
A. Interaction potentials 


B. Reference interaction site model theory 


In the RISM framework^ the pair structure of a fluid 
composed by identical two-site molecules is characterized 
by a set of four site-site intermolecular pair correlation 
functions hij(r) = gij(r) — 1 where (i,j) = (1,2) and 
gij{r) are the site-site radial distribution functions. The 
hij (r) are related to a set of intermolecular direct corre¬ 
lation functions Cij(r) by a matrix generalization of the 
Ornstein-Zernike equation for simple fluidsj^l expressed 
in the /c-space as: 


H(/c) = W{k)C{k)W{k) + pW{k)C{k)H{k) , (2) 


where H = [hij{k)], C = [cij{k)], and W = [wij{k)] 
are 2x2 symmetric matrices; the elements Wij(k) are 
the Fourier transforms of the intramolecular correlation 
functions, written explicitly as: 


{k) 


sm[kLij] 

kLij 


( 3 ) 


where the bond length is given either by Lij = tr, 
if * ^ j, or by Lij = 0 , otherwise. In order to solve 
the RISM equation (I 2 |), we have adopted in this work 
the Kovalenko-Hirata (KH) closure^i^ that assumes — 
beside the exact expression gij(r) = 0 if r < cr — the 
direct correlation functions outside the hard-core to be 
approximated by a combination of the Mean Spherical 
Approximation (MSA) and HyperNetted Chain (HNC) 
expressions^ 


HNC if g^J{r)< I 
MSA if gij (r) > 1 

exp[-/ 3 Vy (r) -k -f^jir)] - -pj{r) - 1 


The sequence of models investigated in this work is 
schematically represented in Fig. [ 1 ] they are constituted 


where /3 = 1 /T* and 7^ (r) = hij{r) — Cij{r). We have 
implemented the numerical solution of the RISM/KH 



3 


scheme by means of a standard iterative Picard algo¬ 
rithm, on a grid of 8192 points with a mesh Ar = 0.005(t. 

In order to calculate the thermodynamic properties, 
we have used standard thermodynamic integration for- 
mulaa^i that already proven to be successful for SW 
dumbbellsi^ Specifically, we have calculated the ex¬ 
cess free energy per particle via the energy route along 
constant-density paths according to: 

/3i^-(/3) /3F-(/3 = 0) , 


the reader to the original works and our recent applica¬ 
tion^ for full details. 


C. Optimized perturbation theory 

The starting point of OPT is given by the optimized 
equation of state recently proposed in Ref. HI, where the 
free energy of a fluid composed by SW particles with well 
depth e and width A, at temperature T and density p is 
expressed as a perturbation series in the form: 


where U{fi) is the internal energy of the fluid at a given 
inverse temperature /3 and = 0 ) corresponds to 

the excess free energy of tangent hard-sphere dumbbells. 
For this, the most classical expressions known in the lit¬ 
erature are an analytic expression due to Tildesley and 
Street^ (and devised to fit their Monte Carlo data) 
and the Wertheim equation of state for tangent hard 
dumbbells;^ reading: 


= 0 ) Up - ipU 

- = n- 


N 


( 1 -? 7)2 


+ {n-l) 


3 ln(l — 77) -I- In 


2-77 


-I- In (677/71) — 1 , 


( 6 ) 


where 77 = 2 ( 7 r/ 6 )p* is the dumbbell packing fraction. 
We have explicitly verified that such equation agrees with 
the free energy expression for hard dumbbells provided 
in Ref. [H. Equation ([6]) has the additional advantage to 
avoid free parameters and the fit of simulation data, so 
that we will use it through all the calculations. 

The excess pressure is then obtained through the (nu¬ 
merical) derivative of the l.h.s. of Eq. (I5|) with respect to 
the density, taken at constant temperature, while coex¬ 
istence conditions are determined by searching for equal 
chemical potentials and pressures (at a given tempera¬ 
ture) of the two coexisting phases. Finally, critical tem¬ 
peratures and densities are obtained by fitting the coexis¬ 
tence points through the scaling law for the densities and 
the law of rectilinear diameters with an effective critical 
exponent (5=0.32.— We refer the reader to our previous 
worb^S for full details on the above procedure. 

Using the symmetric model as a test case, we have also 
investigated the quality of predictions obtained by differ¬ 
ent approaches within the RISM framework. Specifically, 
we have studied the HNC closure, that assumes expres¬ 
sions in the top lines of Eq. o holding for all r > ct 
independently on the value of g{r). Moreover, within 
the KH closure, we have studied a different route from 
structure to thermodynamics, based on closed formulae 
(i.e requiring no thermodynamic integration) to calculate 
free energy and pressure. Such expressions, derived by 
Kovalenko and Hirata in Refs. [13 and [13, closely follow 
the corresponding relations deduced long ago^ within 
the Ornstein-Zernike/HNC scheme for simple fluids and 
later generalized to the RISM/HNC scheme— We refer 


m/3) 

N 


/3F{(3 = 0) ^ 

N ^ ^ 

m—1 



fm[pA) ■ 


( 7 ) 


In Eq. (O F^^{/3 = 0) is the free energy of the reference 
hard-sphere fluid plus the ideal contribution, and fm are 
77 i-th order perturbation terms. In Ref. [13 the series is 
truncated, and fm are calculated up to the fourth-order. 
The theory includes the Barker-Henderson^ result as a 
truncation of the expansion to second order. The third 
and fourth orders, however, are not included as explicit 
higher order corrections to the Barker-Henderson the¬ 
ory, that would require uncontrolled truncations of the 
hierarchy of correlation functions, but as an optimized 
phenomenological calculation of coefficients, hinging on 
simulations results— It turns out that this theory for 
the SW fluid is accurate over nearly the whole density 
range and subcritical (low) temperatures. We refer the 
reader to the original worb^ for detailed expressions of 
fm {m = 1,..., 4), as well as all other necessary machin¬ 
ery. 

The OPT generalization at issue, pertinent to a fluid 
composed by SW dumbbells, is written explicitly as a 
sum of SW contributions, each taking separately into ac¬ 
count the four {i,j = 1,2) interactions of Eq. (P), in the 
form>2^ 


I3F{(3) 

N 


l3Fi0) 

N 


2 





kBT 


fm{p, •^) 


• ( 8 ) 


As visible, each term in the free energy expansion is 
weighted by a temperature reduced by the correspond¬ 
ing well depth Sij, p is the density of sites 1 or 2 , and 
A takes the fixed value 0.5. Given the spherical symme¬ 
try of each site-site interaction, the molecular anisotropy 
is taken into account by the reference free energy, pro¬ 
vided in our case by the Wertheim equation of state for 
tangent hard dumbbells of Eq. This choice gives 

better results for non-convex particles in comparison with 
the possibility to define an effective hard-sphere diameter 
with the corresponding scaling of density.—; it turns to 
be accurate for a moderate number of beads per chain. 

Given Eq. ®, one can obtain the thermodynamic 
properties from usual thermodynamic relations. For in¬ 
stance, the compressibility factor 


Z = 


P 


op 


/3F'^^ 

N 


T 


( 9 ) 
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can also be expressed as a high-temperature expansion, 


Z = 1 + Z-(/3 = 0)+Y.\J2 


i,j—l *- 


knT 


A) 


( 10 ) 

where Zm = rjd/dr] [fm{p, A)]^. 

As in the RISM framework, we have determined the 
gas-liquid phase boundaries by the condition of thermal, 
mechanical, and chemical equilibrium. Numerical meth¬ 
ods have been used to this task, based in this case on a 
combination of the Levenberg-Marquardt algorithm with 
Gauss and steepest-descent methodsi^ 


III. RESULTS 

In Fig. [2] we compare the energy-route KH and OPT 
predictions for symmetric colloidal dumbbells, i.e. with 
£ii = £22 in Eq. O- Results for the pressure (a) and free 
energy (b) are reported along different isotherms across 
the critical temperature T* w 1.60 (estimated in our 
previous work).— KH and OPT predict the same trends 
both for the pressure and free energy; in particular, they 
provide evidence for a supercritical fluid at T* = 1.90, 
whereas a van der Waals loop — heralding a gas-liquid 
phase separation — develops at T* = 1.60 and becomes 
well defined at T* = 1.40. The quantitative agreement 
between the two theories slightly worsens at low tempera¬ 
tures, where OPT curves fall progressively below the cor¬ 
responding KH ones; this circumstance implies a higher 
OPT critical temperature and, correspondingly, a wider 
gas-liquid coexistence curve. In Fig. [2}: we compare KH 
gas-liquid coexistence curves — as obtained by the en¬ 
ergy route and the closed formulm — with Monte Carlo 
(MC) data by usi^ and other authors— OPT predic¬ 
tions, already provided in Ref. are also reported for 
comparison. One can notice that KH reproduces reason¬ 
ably well all features, whereas OPT, as observed, overes¬ 
timates the two-phase region. Such a discrepancy could 
be possibly ascribed to the truncation involved in the 
OPT procedure (see Eq. ([5|)). In this case improvements 
could be achieved by incorporating higher order terms in 
the free energy expansion, so to take into account more 
accurately density fluctuations, along the lines described 
for instance in Ref. [bj. In general, the OPT scheme for 
SW dumbbells shares the same level of accuracy of the 
original approach for the SW atomic fluid)^ the latter 
performing slightly better in the critical region. 

As for HNC, we have verified that at high temper¬ 
atures no significant differences arise with the pressure 
calculated within the KH scheme (see the curves cor¬ 
responding to T* = 1.90 in the inset of Fig. [2h)- On 
the other hand, upon lowering the temperature, the nu¬ 
merical convergence of HNC becomes progressively more 
problematic, thus precluding from further exploration a 
progressively larger interval of T—p conditions. As a con¬ 
sequence, the conditions for gas-liquid coexistence turn 
to be completely unattainable. This issue is even more 





FIG. 2. Pressure (a) and free energy (b) of symmetric col¬ 
loidal dumbbells along different isotherms (see legends): KH 
energy route (lines), KH closed formulae (triangles) and OPT 
(small squares). In the inset, KH and HNC energy-route pres¬ 
sures at T* = 1.90. Panel c: gas-liquid coexistence by KH 
energy-route (KH-ER, squares), KH closed formulae (KH-CF, 
diamonds), OPT (dashed line)^ and MC (circles)Full 
lines are best-fits of KH and simulation points (see text) with 
corresponding critical points (crosses). 
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FIG. 3. KH (lines) and OPT (small squares) pressures of 
colloidal dumbbells with en = 0.3e22. 


pronounced as the molecular anisotropy increases, leav¬ 
ing KH and OPT schemes as the only reliable approaches 
to study the phase separation properties of our models. 

As for KH closed formula, pressures and free energies 
at different temperatures are also reported in Fig. [21 As 
visible, a van der Waals loop in the pressure is observed 
in the same range of temperatures of KH energy route 
and OPT, even if systematically higher and shallower 
estimates are obtained. The same discrepancy appears to 
affect the predictions for the free energy. As a result (see 
Fig. 121:), the liquid branch of phase coexistence appears 
too low in comparison with simulations, KH energy route 
and OPT, whereas the critical temperature and the whole 
gas branch turn to be quite well predicted. Overall, the 
KH energy route appears more accurate and therefore we 
shall adopt such a scheme — together with the differently 
based OPT approach — to investigate the phase behavior 
of other cases with en < £ 22 - 

As visible from Fig. [31 the agreement between KH 
and OPT predictions for the pressure essentially persists 
as £11 is reduced down to 0.3£22. In Ref. [Ij we have 
documented the progressive shrinkage of the coexistence 
curves toward low temperatures as £11 is reduced. As vis¬ 
ible from Fig. 01 where coexistence curves for £11 = 0.7, 
0.5, and 0.3£22 are reported, this feature is reproduced 
by both theories, essentially at the same level of accu¬ 
racy observed for the symmetric case: KH still provides 
reasonable estimates of the gas-liquid coexistence curve, 
with critical temperatures and densities fairly well pre¬ 
dicted, whereas OPT overestimates the amplitude of the 
two-phase region and, accordingly, the critical tempera¬ 
ture. 

As documented in our recent simulation study^i^ in¬ 
creasing the asymmetry of total interaction promotes 
the development of aggregates in the fluid; specifically 




FIG. 4. Gas-liquid coexistence curves of colloidal dumbbells 
upon progressively decreasing en. 
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FIG. 5. Gas-liquid coexistence of colloidal dumbbells with 

£ll = 0.l£22. 


at Ell = 0 . 1 e 22 , dumbbells self-assemble into micelles 
at low densities and planar structures (lamellae) at in¬ 
termediate/high densities, provided the temperature is 
low enough. Moreover the self-assembly process tends to 
destabilize the gas-liquid phase separation - completely 
disappearing for Janus dumbbells, i.e. with eu =0 — 
with the coexistence curve stretched toward very low gas 
densities and a critical temperature drastically lower than 
in the symmetric case. At variance with previous cases, 
both KH and OPT theories fail to reproduce such a be¬ 
havior, largely overestimating the coexistence region at 
Ell = 0.1 e 22, as visible in Fig. [SJ 

Indeed, we may surmise that liquid state theories, ex¬ 
plicitly tailored to describe a homogeneous fluid envi¬ 
ronment, progressively worsen as local inhomogeneities 
tend to develop in the system. Also, the progressive dis¬ 
placement of binodal curves toward lower temperatures 
may negatively affect the OPT predictions based, as de¬ 
scribed, on a high-temperature expansion of the free en¬ 
ergy. 

Numerical values of KH, OPT and MC critical param¬ 
eters for all Ell investigated in this work are collectively 
reported in Tab. H] Critical temperatures and densities 
as functions of Eu are shown in Fig. [6] Remarkably, 
both theories agree with simulations in providing a lin¬ 
ear dependence of the critical temperature on Eh. Such a 
circumstance is compatible with a mean held behavior^ 
with the only exception (for MC data) of the Eu = O.IE 22 
case. Conversely, the critical density p* follows quite a 
different trend, keeping an almost constant value from 
the symmetric case down to Eu = 0.5 e 22, with a marked 
drop upon further lowering this parameter. Surprisingly 
enough, OPT qualitatively reproduces such a behavior, 
actually predicting a critical density close to simulation 
data for Eh = O.IE 22 . On the other hand, the dependence 


TABLE I. Critical temperatures, densities and chemical po¬ 
tentials for colloidal dumbbells with variable £11 as obtained 
by MC,-i^ KH and OPT. 


£ii/e22 


T* 

Pc 

Mc/£22 

1.0 

MC 

1.566 

0.146 

-5.680 


KH 

1.616 

0.135 

-5.973 


OPT 

1.842 

0.106 

-6.549 

0.7 

MC 

1.213 

0.146 

-4.524 


KH 

1.252 

0.133 

-4.610 


OPT 

1.430 

0.105 

-5.025 

0.5 

MC 

0.980 

0.143 

-3.725 


KH 

1.017 

0.130 

-3.762 


OPT 

1.162 

0.103 

-4.107 

0.3 

MC 

0.739 

0.125 

-2.997 


KH 

0.794 

0.126 

-2.964 


OPT 

0.905 

0.097 

-3.248 

0.1 

MC 

0.466 

0.068 

-2.635 


KH 

0.591 

0.116 

-2.284 


OPT 

0.682 

0.080 

-2.205 


of KH p* on Ell is less pronounced, failing to reproduce 
in particular the very low values attained at Eu = O.IE 22 . 

In a recent study,— we have analyzed the RISM per¬ 
formances for a close class of SW dumbbells, character¬ 
ized for instance by an attraction range A = 0.1, shorter 
than that investigated here (namely A = 0.5). In that 
case we have shown that the HNC theory is able to de¬ 
scribe the low temperature regime whereupon the gas- 
liquid coexistence takes place; in particular, the theory 
predicts the disappearance of a stable phase separation 
for Ell < 0.4 e 22, i.e. at a value slightly larger than that 
observed in simulations (eh < O. 2 E 22 ). Also for A = 0.1 
RISM agrees with simulations in predicting the linear 
scaling of the critical temperature discussed in Fig. [6^, 
even if with a larger discrepancy with respect to the 
present case; a better agreement is instead observed for 
the critical density, essentially because this quantity, at 
A = 0.1, keeps an almost constant value as a function 
of Ell, at variance with the sudden drop documented in 

Fig.ini). 

IV. CONCLUSIONS 

We have presented a theoretical study of thermody¬ 
namic properties and gas-liquid phase coexistence of 
anisotropic colloidal dumbbells, where the anisotropy 
stems from the fact that, although the two colloids form¬ 
ing the dumbbell are identical in sizes, their interactions 
are not, as one of the two attractive square-well depth 
is gradually reduced to zero toward the Janus dumbbells 
limit. We have discussed and contrasted to each other 
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FIG. 6. Critical temperatures (a) and densities (b) of colloidal 
dumbbells as functions of en. 


the performances of two venerable and well established 
tools, the RISM and OPT theories of liquids. As for the 
RISM, this theory hinges on the solution of the integral 
equations of liquids, thus providing the intermolecular 
correlations, whereby all thermodynamic properties can 
be obtained. Using the symmetric model as a test case, 
we have shown that the Kovalenko-Hirata closure (com¬ 
bining efficiently MSA and HNC approximations) com¬ 
plemented with the energy route from structure to ther¬ 
modynamics is superior to the HNC approximation and 
other procedures to calculate the thermodynamic prop¬ 
erties. As for the OPT theory, this method hinges on a 
fourth-order high-temperature perturbative expansion of 
the free energy of the square-well fluid, recently general¬ 
ized to deal with molecular fluids. 

The accuracy of both methods has been assessed by 
a direct comparison with recent numerical simulations. 
We have found that both approaches provide reasonable 


predictions, with RISM slightly outperforming OPT in 
almost all cases with the notable exception of the critical 
density, where OPT is more accurate. Latter predictions 
could be improved by finding a way to fully include in 
the formalism the intramolecular correlations, presently 
taken into account only at the level of the reference free 
energy. 

We emphasize the importance of having a reliable the¬ 
oretical tool to provide a quick estimate of thermody¬ 
namics of colloidal dumbbells. Recent experiments;^ as 
well as complementary numerical work;i^ show that in 
most of the interesting cases the two colloids forming the 
dumbbells are asymmetric both in size and interaction. 
This means that a huge number of possible combinations 
could be in principle studied to probe many topologically 
different resulting phase diagrams. A reliable theoretical 
approach would then be invaluable in pinning down the 
most interesting cases where a more sophisticated, and 
computationally more demanding, calculation could be 
performed. 
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